---
title: "Trust after Tragedy: How Traumatic Events Impact Trust in Government(s)"
author: "Wayde Z.C. Marsh"
date: "07/28/25"
output:
  html_document:
    toc: true
    code_folding: hide
  latex_engine: xelatex
  number_sections: yes
mainfont: Helvetica
---

# Setting up the Data {.tabset}

## Loading packages and data

```{r setup, message = FALSE, warning = FALSE, results ='hide'}

##### Packages #####

library(descr)
library(ggplot2)
library(reshape2)
library(psych)
library(GPArotation)
library(dplyr)
library(foreign)
library(ggpubr)
library(ggeffects)
library(data.table)
library(haven)
library(lavaan)
library(tidyverse)
library(modelsummary)
library(fixest)
library(tidySEM)

setwd('[file path here]')
knitr::opts_knit$set(root.dir = '[file path here]')

```

# Load Data

```{r}

load('study1.Rda') # study I
load('study2.Rda') # study II
load('study3.Rda') # study III
load('study4_corrected.Rda') # study IV
load('study5.Rda') # study V

dat2$study <- rep('Two', nrow(dat2))
dat3$study <- rep('Three', nrow(dat3))
dat4$study <- rep('Four', nrow(dat4))
dat5$study <- rep('Five', nrow(dat5))


datp <- do.call('rbind', list(dat2, dat3, dat4, dat5))

```

# Factor Analysis

```{r FA}

fa_vars <- c('ptsi_1', 'ptsi_2', 'ptsi_3', 'ptsi_4', 'ptsi_5',
             'ptgi_1', 'ptgi_2', 'ptgi_3', 'ptgi_4', 'ptgi_5')

fa1 <- dat1[fa_vars]
fa2 <- dat2[fa_vars]
fa3 <- dat3[fa_vars]
fa4 <- dat4[fa_vars]
fa5 <- dat5[fa_vars]

eigen(cor(na.omit(fa1))) # two eigen values > 2 for all confirming 2 factors
eigen(cor(na.omit(fa2)))
eigen(cor(na.omit(fa3)))
eigen(cor(na.omit(fa4)))
eigen(cor(na.omit(fa5)))

## Parallel Analysis

sem.parallel<-function(x, R=1000, percent=.95){
  N<-nrow(x)
  P<-ncol(x)
  
  # original eigenvalues
  orig.eigen<-eigen(cor(na.omit(x)))$values
  
  pa.eigen<-NULL
  
  for (i in 1:R){
    temp<-array(rnorm(N*P), c(N,P))
    temp.eigen<-eigen(cor(temp))
    pa.eigen<-rbind(pa.eigen, temp.eigen$values)
  }
  
  m.eigen<-apply(pa.eigen, 2, mean)
  p.eigen<-apply(pa.eigen, 2, quantile, prob=percent)
  ## plot the eigen values
  plot(orig.eigen, type='b', xlab='Number of factors', ylab='Eigenvalues')
  lines(m.eigen, lty=3)
  lines(p.eigen, lty=2)
  abline(h = 1, col = 'red')
  legend('topright', c('Data', 'Average', paste(percent*100, '%', sep='')), lty=c(1, 3, 2))
  invisible(list(eigen=orig.eigen, pa.eigen=pa.eigen))
}

set.seed(52617)

sem.parallel(fa1) # Parallel Analysis confirms 2 factors for all dfs
sem.parallel(fa2)
sem.parallel(fa3)
sem.parallel(fa4)
sem.parallel(fa5)

# EFA

### EFA ####
fit1 <- fa(fa1, nfactors = 2, rotate = "oblimin", fm = "minres", alpha = 0.05)
fit2 <- fa(fa2, nfactors = 2, rotate = "oblimin", fm = "minres", alpha = 0.05)
fit3 <- fa(fa3, nfactors = 2, rotate = "oblimin", fm = "minres", alpha = 0.05)
fit4 <- fa(fa4, nfactors = 2, rotate = "oblimin", fm = "minres", alpha = 0.05)
fit5 <- fa(fa5, nfactors = 2, rotate = "oblimin", fm = "minres", alpha = 0.05)

fit1s <- fa(fa1[, 1:5], nfactors = 2, rotate = "oblimin", fm = "minres", alpha = 0.05)
fit2s <- fa(fa2[, 1:5], nfactors = 2, rotate = "oblimin", fm = "minres", alpha = 0.05)
fit3s <- fa(fa3[, 1:5], nfactors = 2, rotate = "oblimin", fm = "minres", alpha = 0.05)
fit4s <- fa(fa4[, 1:5], nfactors = 2, rotate = "oblimin", fm = "minres", alpha = 0.05)
fit5s <- fa(fa5[, 1:5], nfactors = 2, rotate = "oblimin", fm = "minres", alpha = 0.05)

fit1g <- fa(fa1[, 6:10], nfactors = 2, rotate = "oblimin", fm = "minres", alpha = 0.05)
fit2g <- fa(fa2[, 6:10], nfactors = 2, rotate = "oblimin", fm = "minres", alpha = 0.05)
fit3g <- fa(fa3[, 6:10], nfactors = 2, rotate = "oblimin", fm = "minres", alpha = 0.05)
fit4g <- fa(fa4[, 6:10], nfactors = 2, rotate = "oblimin", fm = "minres", alpha = 0.05)
fit5g <- fa(fa5[, 6:10], nfactors = 2, rotate = "oblimin", fm = "minres", alpha = 0.05)

print(fit1$loadings, cutoff = 0.6)
print(fit2$loadings, cutoff = 0.6)
print(fit3$loadings, cutoff = 0.6)
print(fit4$loadings, cutoff = 0.6)
print(fit5$loadings, cutoff = 0.6)

psych::alpha(fa1[, 1:5])
psych::alpha(fa2[, 1:5])
psych::alpha(fa3[, 1:5])
psych::alpha(fa4[, 1:5])
psych::alpha(fa5[, 1:5])

psych::alpha(fa1[, 6:10])
psych::alpha(fa2[, 6:10])
psych::alpha(fa3[, 6:10])
psych::alpha(fa4[, 6:10])
psych::alpha(fa5[, 6:10])

```

# Descriptives

```{r Descriptives}

datp$tr_type <- ifelse(datp$trauma1_1 == 1, 'Sexual Violence',
                ifelse(datp$trauma1_2 == 1, 'Natural Disaster',
                ifelse(datp$trauma1_3 == 1, 'Industrial Disaster',
                ifelse(datp$trauma1_4 == 1, 'School Shooting',
                ifelse(datp$trauma1_5 == 1, 'Church Shooting',
                ifelse(datp$trauma1_6 == 1, 'Other Mass Shooting',
                ifelse(datp$trauma1_7 == 1, 'Car Accident',
                ifelse(datp$trauma1_8 == 1, 'Plane Accident',
                ifelse(datp$trauma1_9 == 1, 'Assault',
                ifelse(datp$trauma1_10 == 1, 'Warfare',
                ifelse(datp$trauma1_11 == 1, 'Other', 'None')))))))))))

datp$rtr_type <- ifelse(datp$trauma2_1 == 1, 'Sexual Violence',
                 ifelse(datp$trauma2_2 == 1, 'Natural Disaster',
                 ifelse(datp$trauma2_3 == 1, 'Industrial Disaster',
                 ifelse(datp$trauma2_4 == 1, 'School Shooting',
                 ifelse(datp$trauma2_5 == 1, 'Church Shooting',
                 ifelse(datp$trauma2_6 == 1, 'Other Mass Shooting',
                 ifelse(datp$trauma2_7 == 1, 'Car Accident',
                 ifelse(datp$trauma2_8 == 1, 'Plane Accident',
                 ifelse(datp$trauma2_9 == 1, 'Assault',
                 ifelse(datp$trauma2_10 == 1, 'Warfare',
                 ifelse(datp$trauma2_11 == 1, 'Other', 'None')))))))))))

datp <- datp %>%
  mutate(tr_type_simp = recode(tr_type, 
                               'Sexual Violence' = 4,
                               'Natural Disaster' = 1,
                               'Industrial Disaster' = 2,
                               'School Shooting' = 4,
                               'Church Shooting' = 4,
                               'Other Mass Shooting' = 4,
                               'Car Accident' = 3,
                               'Plane Accident' = 3,
                               'Assault' = 4,
                               'Warfare' = 4,
                               'Other' = 2,
                               'None' = 0))

datp <- datp %>%
  mutate(rtr_type_simp = recode(rtr_type, 
                               'Sexual Violence' = 4,
                               'Natural Disaster' = 1,
                               'Industrial Disaster' = 2,
                               'School Shooting' = 4,
                               'Church Shooting' = 4,
                               'Other Mass Shooting' = 4,
                               'Car Accident' = 3,
                               'Plane Accident' = 3,
                               'Assault' = 4,
                               'Warfare' = 4,
                               'Other' = 2,
                               'None' = 0))


```

## Figures

```{r}

plot_type_df <- data.frame(type = c('aNone', 'bNatural Disaster', 'cIndustrial Disaster',
                                    'dCar Accident', 'ePlane Accident', 
                                    'fSchool Shooting', 'gChurch Shooting',
                                    'hOther Mass Shooting', 'iAssault', 'jWarfare',
                                    'kSexual Violence', 'lOther'),
                           experienced = c(rep('Respondent', 12), 
                                           rep('Friend/Family Member', 12)),
                           exposed = c(41.0, 21.4, 1.5, 8.4, 0.2, 0.8, 0.2, 0.4, 1.6, 0.9,
                                    21.7, 2.0, 36.6, 16.6, 1.1, 8.1, 0.2, 0.4, 0.1, 0.5, 2.1,
                                    2.9, 30.2, 1.0))

p_type1 <- ggplot(plot_type_df, aes(x = type, y = exposed, fill = experienced)) +
  theme_bw() +
  geom_bar(stat = 'identity', position = position_dodge()) +
  scale_fill_manual(name = 'Who Experienced',
                    values = c('grey65', 'gray30')) +
  scale_x_discrete(name = 'Type of Traumatic Event', 
                   labels = c('None', 'Natural Disaster', 'Industrial Disaster',
                                    'Car Accident', 'Plane Accident', 
                                    'School Shooting', 'Church Shooting',
                                    'Other Mass Shooting', 'Assault', 'Warfare',
                                    'Sexual Violence', 'Other')) +
  ylim(c(0, 50)) +
  ylab('Percentage of Sample Experienced') +
  labs(title = 'Sexual Violence and Natural Disasters more Common Traumatic Events Experienced') +
  theme(axis.text=element_text(size = 10)) + 
  theme(plot.title = element_text(hjust = 0.5)) + 
    theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1)) +
  theme(plot.caption = element_text(hjust = 0)) +
  theme(text=element_text(size = 10, family="serif"))

# ggsave('p_type1.png', p_type1, 'png', '~/Dropbox/Wayde/Dissertation/Trust',
#     width = 190, height = 100, units = 'mm')

plot_type2_df <- data.frame(type = c('aNone', 'bNatural Disaster', 'cAnthropogenic Disaster',
                                    'dAnthropogenic Interpersonal Accidental', 
                                    'eAnthropogenic Interpersonal Intentional'),
                           experienced = c(rep('Respondent', 5), 
                                           rep('Friend/Family Member', 5)),
                           exposed = c(41.0, 21.4, 3.5, 8.5, 25.5,
                                       36.6, 16.6, 2.1, 8.3, 36.3))

p_type2 <- ggplot(plot_type2_df, aes(x = type, y = exposed, fill = experienced)) +
  theme_bw() +
  geom_bar(stat = 'identity', position = position_dodge()) +
  scale_fill_manual(name = 'Who Experienced',
                    values = c('grey65', 'gray30')) +
  scale_x_discrete(name = '\nType of Traumatic Event', 
                   labels = c('None', 'Natural Disaster', 'Anthropogenic \nDisaster',
                                    'Anthropogenic \nInterpersonal \nAccidental', 
                                    'Anthropogenic \nInterpersonal \nIntentional')) +
  ylim(c(0, 50)) +
  ylab('Percentage of Sample Experienced') +
  labs(title = 'Prevalence of Exposure to Traumatic Events, by Type') +
  theme(axis.text=element_text(size = 10)) + 
  theme(plot.title = element_text(hjust = 0.5)) + 
    theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1)) +
  theme(plot.caption = element_text(hjust = 0)) +
  theme(text=element_text(size = 10, family="serif"))

# ggsave('p_type2.png', p_type2, 'png', '~/Dropbox/Wayde/Dissertation/Trust',
#     width = 190, height = 100, units = 'mm')

## PTS and PTG by type

df_tr <- datp %>%
  group_by(tr_type_simp) %>%
  summarize(avg = mean(ptsi, na.rm = T))

df_tr2 <- datp %>%
  group_by(tr_type_simp) %>%
  summarize(avg = mean(ptgi, na.rm = T))

df_rtr <- datp %>%
  group_by(rtr_type_simp) %>%
  summarize(avg = mean(ptsi, na.rm = T))

df_rtr2 <- datp %>%
  group_by(rtr_type_simp) %>%
  summarize(avg = mean(ptgi, na.rm = T))

df_tr <- rbind(df_tr, df_tr2)
df_rtr <- rbind(df_rtr, df_rtr2)

df_tr$exp <- rep('Respondent', 10)
df_rtr$exp <- rep('Close Friend/Family Member', 10)

df_tr <- subset(df_tr, tr_type_simp > 0)
df_rtr <- subset(df_rtr, rtr_type_simp > 0)

df_tr$psych <- c(rep('PTS', 4), rep('PTG', 4))
df_rtr$psych <- c(rep('PTS', 4), rep('PTG', 4))

names(df_rtr) <- names(df_tr)

df_tr <- rbind(df_tr, df_rtr)

p_ptsptg <- ggplot(df_tr, aes(x = as.factor(tr_type_simp), y = avg, fill = psych)) +
  theme_bw() +
  facet_wrap(vars(exp)) +
  geom_bar(stat = 'identity', position = position_dodge()) +
  scale_fill_manual(name = 'Psychological \nResponse',
                    values = c('grey65', 'gray30')) +
  scale_x_discrete(name = '\nType of Traumatic Event', 
                   labels = c('Natural Disaster', 'Anthropogenic \nDisaster',
                                    'Anthropogenic \nInterpersonal \nAccidental', 
                                    'Anthropogenic \nInterpersonal \nIntentional')) +
  ylim(c(0, 1)) +
  ylab('Mean Posttraumatic Response') +
  # labs(title = 'Prevalence of Exposure to Traumatic Events, by Type') +
  theme(axis.text=element_text(size = 10)) + 
  theme(plot.title = element_text(hjust = 0.5)) + 
    theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1)) +
  theme(plot.caption = element_text(hjust = 0)) +
  theme(text=element_text(size = 10, family="serif"))

# ggsave('p_ptsptg.png', p_ptsptg, 'png', '~/Dropbox/Wayde/Dissertation/Trust',
#     width = 190, height = 100, units = 'mm')

df_bl <- datp %>%
  group_by(black) %>%
  summarize(avg = mean(ptsi, na.rm = T))

df_bl2 <- datp %>%
  group_by(black) %>%
  summarize(avg = mean(ptgi, na.rm = T))

df_hi <- datp %>%
  group_by(hisp) %>%
  summarize(avg = mean(ptsi, na.rm = T))

df_hi2 <- datp %>%
  group_by(hisp) %>%
  summarize(avg = mean(ptgi, na.rm = T))

df_fem <- datp %>%
  group_by(female) %>%
  summarize(avg = mean(ptsi, na.rm = T))

df_fem2 <- datp %>%
  group_by(female) %>%
  summarize(avg = mean(ptgi, na.rm = T))

df_bl <- data.frame(
  psych = c(rep('PTS', 2), rep('PTG', 2)),
  Race = c('Non-Black', 'Black', 'Non-Black', 'Black'),
  avg = c(0.300, 0.375, 0.497, 0.658)
)

df_hi <- data.frame(
  psych = c(rep('PTS', 2), rep('PTG', 2)),
  Race = c('Non-Hispanic', 'Hispanic', 'Non-Hispanic', 'Hispanic'),
  avg = c(0.302, 0.386, 0.503, 0.588)
)

df_fem <- data.frame(
  psych = c(rep('PTS', 2), rep('PTG', 2)),
  Race = c('Non-Female', 'Female', 'Non-Female', 'Female'),
  avg = c(0.258, 0.347, 0.468, 0.550)
)

df_demo <- do.call('rbind', list(df_bl, df_hi, df_fem))

df_demo$demo <- c(rep('Race', 4), rep('Ethnicity', 4), rep('Sex', 4))

p_demo <- ggplot(df_demo, aes(x = as.factor(psych), y = avg, fill = Race)) +
  theme_bw() +
  facet_wrap(vars(demo)) +
  geom_bar(stat = 'identity', position = position_dodge()) +
  scale_fill_manual(name = 'Demographic',
                    values = c('gray5', 'gray5', 'gray5',
                               'grey65', 'grey65', 'grey65')) +
  # scale_x_discrete(name = '\nType of Traumatic Event',
  #                  labels = c('Natural Disaster', 'Anthropogenic \nDisaster',
  #                                   'Anthropogenic \nInterpersonal \nAccidental',
  #                                   'Anthropogenic \nInterpersonal \nIntentional')) +
  xlab('Psychological Response') +
  ylim(c(0, 1)) +
  ylab('Mean Posttraumatic Response') +
  # labs(title = 'Prevalence of Exposure to Traumatic Events, by Type') +
  theme(axis.text=element_text(size = 10)) + 
  theme(plot.title = element_text(hjust = 0.5)) + 
    theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust=1)) +
  theme(plot.caption = element_text(hjust = 0)) +
  theme(text=element_text(size = 10, family="serif"))

# ggsave('p_demo.png', p_demo, 'png', '~/Dropbox/Wayde/Dissertation/Trust',
#     width = 190, height = 100, units = 'mm')

```

## Tests

### Regression

```{r}

# demographics
## income less than 0.3 is lower class, upper class is greater than 0.8 (YouGov scale)
## corresponds to less than $50k a year and more than $120k a year (family income)
datp$class <- ifelse(datp$income < 0.3, 0,
                ifelse(datp$income > 0.8, 1, 0.5))
datp$class2 <- ifelse(datp$income < 0.3, 'lower',
                ifelse(datp$income > 0.8, 'upper', 'middle'))
dat1$class <- ifelse(dat1$income < 0.3, 0,
                ifelse(dat1$income > 0.8, 1, 0.5))
dat1$class2 <- ifelse(dat1$income < 0.3, 'lower',
                ifelse(dat1$income > 0.8, 'upper', 'middle'))
dat1$lower <- ifelse(dat1$class == 0, 1, 0)

datp2 <- subset(datp, tr_type_simp > 0 | rtr_type_simp > 0)
datp2$lower <- ifelse(datp2$class == 0, 1, 0)
datp$lower <- ifelse(datp$class == 0, 1, 0)

datp$Direct <- ifelse(datp$tr_exp == 1, 1, 0)
datp$Indirect <- ifelse(datp$rtr_exp == 1, 1, 0)

datp$tr_exp_t <-ifelse(is.na(datp$tr_exp_t), 50, datp$tr_exp_t)
datp$rtr_exp_t <-ifelse(is.na(datp$rtr_exp_t), 50, datp$rtr_exp_t)

datp$Time <- ifelse(datp$tr_exp == 1 & datp$rtr_exp != 1, datp$tr_exp_t,
             ifelse(datp$rtr_exp == 1 & datp$tr_exp != 1, datp$rtr_exp_t,
             ifelse(datp$tr_exp == 1 & datp$rtr_exp == 1, pmin(datp$tr_exp_t, datp$rtr_exp_t), NA)))

datp$Time2 <- ifelse(datp$Time == 0, 1,
              ifelse(datp$Time == 1, 4,
              ifelse(datp$Time < 0.5 & datp$Time > 0, 2,
              ifelse(datp$Time > 0.5 & datp$Time < 1, 3, NA))))

# datp$Time <- ifelse(datp$Time == 50 | datp$Time == 50, NA, datp$Time)

datp$PTS <- datp$ptsi
datp$PTG <- datp$ptgi
datp$Black <- datp$black
datp$Female <- datp$female
datp$Class <- datp$class
datp$Republican <- datp$gop
datp$PolInterest <- datp$newsint

datp$rtr_type_simp <- ifelse(is.na(datp$rtr_type_simp), 0, datp$rtr_type_simp)
datp$Type <- ifelse(datp$tr_exp == 1 & datp$rtr_exp != 1, datp$tr_type_simp,
             ifelse(datp$rtr_exp == 1 & datp$tr_exp != 1, datp$rtr_type_simp,
             ifelse((datp$tr_exp == 1 & datp$rtr_exp == 1) & datp$tr_exp_t > datp$rtr_exp_t, datp$rtr_type_simp, 
             ifelse((datp$tr_exp == 1 & datp$rtr_exp == 1) & datp$tr_exp_t <= datp$rtr_exp_t, datp$tr_type_simp, NA))))

# save(datp, file = 'datp_final.Rda')

m_type1 <- feols(ptgi ~ factor(tr_type_simp) + tr_exp_t + black + hisp + female + lower | study,
              weights = ~weight,
              data = subset(datp2, tr_type_simp > 0))

m_type2 <- feols(ptsi ~ factor(tr_type_simp) + tr_exp_t + black + hisp + female + lower | study,
              weights = ~weight,
              data = subset(datp2, tr_type_simp > 0))

m_type1_r <- feols(ptgi ~ factor(rtr_type_simp) + rtr_exp_t + black + hisp +female + lower | study,
              weights = ~weight,
              data = subset(datp2, rtr_type_simp > 0))

m_type2_r <- feols(ptsi ~ factor(rtr_type_simp) + rtr_exp_t + black + hisp + female + lower | study,
              weights = ~weight,
              data = subset(datp2, rtr_type_simp > 0))

summary(m_type1); summary(m_type2); summary(m_type1_r); summary(m_type2_r)

modelsummary(models = list('Direct PTG'  = m_type1,
                           'Direct PTS' = m_type2,
                           'Indirect PTG' = m_type1_r,
                           'Indirect PTS'  = m_type2_r),
             output = 'latex',
             stars = c('*' = .05),
             gof_omit = 'AIC|BIC|F|RMSE|Log.Lik.',
             escape = FALSE)


m_demo1 <- feols(ptsi ~ black + hisp + female + lower,
                 data = dat1)
m_demo2 <- feols(ptgi ~ black + hisp + female + lower,
                 data = dat1)

m_demo3 <- feols(ptsi ~ black + hisp + female + lower | study,
                 weights = ~weight,
                 data = datp2)
m_demo4 <- feols(ptgi ~ black + hisp + female + lower | study,
                 weights = ~weight,
                 data = datp2)
summary(m_demo1); summary(m_demo2)

modelsummary(models = list('Study 1 PTG'  = m_demo2,
                           'Study 1 PTS' = m_demo1,
                           'Studies 2-5 PTG'  = m_demo4,
                           'Studies 2-5 PTS' = m_demo3),
             output = 'latex',
             stars = c('*' = .05),
             gof_omit = 'AIC|BIC|F|RMSE|Log.Lik.',
             escape = FALSE)

# m_trst1 <- feols(trstgov ~ ptsi + black + female + income | study,
#                  weights = ~weight,
#               data = subset(datp, tr_type_simp > 0))
# summary(m_trst1)


## MAIN ANALYSES ####

m_trst <- feols(c(trstgov, trststate, trstloc, trstgen, trstgen2) ~ PTS + PTG + Time +
                  Black + Female + Class + Republican + PolInterest | study,
                weights = ~weight,
                data = datp)
summary(m_trst)

m_trst_direct <- feols(c(trstgov, trststate, trstloc, trstgen, trstgen2) ~ 
                  PTS + PTG + Time +
                  Black + Female + Class + Republican + PolInterest | study,
                weights = ~weight,
                data = subset(datp, Direct == 1))
summary(m_trst_direct)

m_trst_indirect <- feols(c(trstgov, trststate, trstloc, trstgen, trstgen2) ~ 
                  PTS + PTG + Time +
                  Black + Female + Class + Republican + PolInterest | study,
                weights = ~weight,
                data = subset(datp, Indirect == 1))
summary(m_trst_indirect)

m_vote <- feols(c(presvote16, presvote20) ~ PTS + PTG + Time +
                  Black + Female + Class + Republican + PolInterest | study,
                weights = ~weight,
                data = datp)
summary(m_vote)

m_vote_direct <- feols(c(presvote16, presvote20) ~ PTS + PTG + Time +
                  Black + Female + Class + Republican + PolInterest | study,
                weights = ~weight,
                data = subset(datp, Direct == 1))
summary(m_vote_direct)

m_vote_indirect <- feols(c(presvote16, presvote20) ~ PTS + PTG + Time +
                  Black + Female + Class + Republican + PolInterest | study,
                weights = ~weight,
                data = subset(datp, Indirect == 1))
summary(m_vote_indirect)

modelsummary(models = list('Federal'  = m_trst[1],
                           'State'  = m_trst[2],
                           'Local'  = m_trst[3],
                           'Interpersonal I'  = m_trst[4],
                           'Interpersonal II'  = m_trst[5]),
             output = 'latex',
             stars = c('*' = .05),
             gof_omit = 'AIC|BIC|F|RMSE|Log.Lik.',
             escape = FALSE)

modelsummary(models = list('Federal'  = m_trst_direct[1],
                           'Federal'  = m_trst_indirect[1],
                           'State'  = m_trst_direct[2],
                           'State'  = m_trst_indirect[2],
                           'Local'  = m_trst_direct[3],
                           'Local'  = m_trst_indirect[3],
                           'Interpersonal I'  = m_trst_direct[4],
                           'Interpersonal I'  = m_trst_indirect[4],
                           'Interpersonal II'  = m_trst_direct[5],
                           'Interpersonal II'  = m_trst_indirect[5]),
             output = 'latex',
             stars = c('*' = .05),
             gof_omit = 'AIC|BIC|F|RMSE|Log.Lik.',
             escape = FALSE)

modelsummary(models = list('2016 Voting'  = m_vote[1],
                           '2020 Voting'  = m_vote[2]),
             output = 'latex',
             stars = c('*' = .05),
             gof_omit = 'AIC|BIC|F|RMSE|Log.Lik.',
             escape = FALSE)

modelsummary(models = list('2016 Voting'  = m_vote_direct[1],
                           '2016 Voting'  = m_vote_indirect[1],
                           '2020 Voting'  = m_vote_direct[2],
                           '2020 Voting'  = m_vote_indirect[2]),
             output = 'latex',
             stars = c('*' = .05),
             gof_omit = 'AIC|BIC|F|RMSE|Log.Lik.',
             escape = FALSE)

# all
coefplot(m_trst,main = '',
                   col = 1,
                   drop = c('Time', 'Black', 'Female', 'Class', 'Republican', 'PolInterest'),
                   value.lab = 'Estimate with __ci__ Confidence Interval',
                   par(family = 'serif'))
legend('bottomright', xpd = T,
         col = 1, pch = c(20, 17, 15, 21, 24),
       cex = 0.7,
       legend = c('Federal', 'State', 'Local', 
                              'Interpersonal I', 'Interpersonal II'),
       title = 'Trust DV')

# direct
coefplot(m_trst_direct,main = '',
                   col = 1,
                   drop = c('Time', 'Black', 'Female', 'Class', 'Republican', 'PolInterest'),
                   value.lab = 'Estimate with __ci__ Confidence Interval')
legend('bottomright', xpd = T,
         col = 1, pch = c(20, 17, 15, 21, 24),
       cex = 0.7,
       legend = c('Federal', 'State', 'Local', 
                              'Interpersonal I', 'Interpersonal II'),
       title = 'Trust DV')

# indirect
coefplot(m_trst_indirect,main = '',
                   col = 1,
                   drop = c('Time', 'Black', 'Female', 'Class', 'Republican', 'PolInterest'),
                   value.lab = 'Estimate with __ci__ Confidence Interval')
legend('bottomright', xpd = T,
         col = 1, pch = c(20, 17, 15, 21, 24),
       cex = 0.7,
       legend = c('Federal', 'State', 'Local', 
                              'Interpersonal I', 'Interpersonal II'),
       title = 'Trust DV')

# VOTING

#all
coefplot(m_vote,main = '',
                   col = 1,
                   drop = c('Time', 'Black', 'Female', 'Class', 'Republican', 'PolInterest'),
                   sep = 0.5,
                   value.lab = 'Estimate with __ci__ Confidence Interval')
legend('bottomright', 
         col = 1, pch = c(20, 17),
       cex = 0.7,
       legend = c('2016 Voting', '2020 Voting'),
       title = 'Turnout DV')

# direct
coefplot(m_vote_direct,main = '',
                   col = 1,
                   drop = c('Time', 'Black', 'Female', 'Class', 'Republican', 'PolInterest'),
                   sep = 0.5,
                   value.lab = 'Estimate with __ci__ Confidence Interval')
legend('bottomright', 
         col = 1, pch = c(20, 17),
       cex = 0.7,
       legend = c('2016 Voting', '2020 Voting'),
       title = 'Turnout DV')

# indirect
coefplot(m_vote_indirect,main = '',
                   col = 1,
                   drop = c('Time', 'Black', 'Female', 'Class', 'Republican', 'PolInterest'),
                   sep = 0.5,
                   value.lab = 'Estimate with __ci__ Confidence Interval')
legend('bottomright', 
         col = 1, pch = c(20, 17),
       cex = 0.7,
       legend = c('2016 Voting', '2020 Voting'),
       title = 'Turnout DV')

#### Type Analysis ####

type_trust1 <- feols(c(trstgov, trststate, trstloc, trstgen, trstgen2) ~ 
                  PTS + PTG + Time +
                  Black + Female + Class + Republican + PolInterest | study,
                weights = ~weight,
                data = subset(datp, Type == 1))
summary(type_trust1)

type_trust4 <- feols(c(trstgov, trststate, trstloc, trstgen, trstgen2) ~ 
                  PTS + PTG + Time +
                  Black + Female + Class + Republican + PolInterest | study,
                weights = ~weight,
                data = subset(datp, Type == 4))
summary(type_trust4)

type_vote1 <- feols(c(presvote16, presvote20) ~ PTS + PTG + Time +
                  Black + Female + Class + Republican + PolInterest | study,
                weights = ~weight,
                data = subset(datp, Type == 1))
summary(type_vote1)

type_vote4 <- feols(c(presvote16, presvote20) ~ PTS + PTG + Time +
                  Black + Female + Class + Republican + PolInterest | study,
                weights = ~weight,
                data = subset(datp, Type == 4))
summary(type_vote4)

modelsummary(models = list('Federal'  = type_trust1[1],
                           'State'  = type_trust1[2],
                           'Local'  = type_trust1[3],
                           'Interpersonal I'  = type_trust1[4],
                           'Interpersonal II'  = type_trust1[5]),
             output = 'latex',
             stars = c('*' = .05),
             gof_omit = 'AIC|BIC|F|RMSE|Log.Lik.',
             escape = FALSE)

modelsummary(models = list('Federal'  = type_trust4[1],
                           'State'  = type_trust4[2],
                           'Local'  = type_trust4[3],
                           'Interpersonal I'  = type_trust4[4],
                           'Interpersonal II'  = type_trust4[5]),
             output = 'latex',
             stars = c('*' = .05),
             gof_omit = 'AIC|BIC|F|RMSE|Log.Lik.',
             escape = FALSE)

modelsummary(models = list('2016 Voting'  = type_vote1[1],
                           '2020 Voting'  = type_vote1[2],
                           '2016 Voting'  = type_vote4[1],
                           '2020 Voting'  = type_vote4[2]),
             output = 'latex',
             stars = c('*' = .05),
             gof_omit = 'AIC|BIC|F|RMSE|Log.Lik.',
             escape = FALSE)

coefplot(type_trust1,main = '',
                   col = 1,
                   drop = c('Time', 'Black', 'Female', 'Class', 'Republican', 'PolInterest'),
                   value.lab = 'Estimate with __ci__ Confidence Interval')
legend('bottomright', xpd = T,
         col = 1, pch = c(20, 17, 15, 21, 24),
       cex = 0.7,
       legend = c('Federal', 'State', 'Local', 
                              'Interpersonal I', 'Interpersonal II'),
       title = 'Trust DV')

coefplot(type_trust4,main = '',
                   col = 1,
                   drop = c('Time', 'Black', 'Female', 'Class', 'Republican', 'PolInterest'),
                   value.lab = 'Estimate with __ci__ Confidence Interval')
legend('bottomright', xpd = T,
         col = 1, pch = c(20, 17, 15, 21, 24),
       cex = 0.7,
       legend = c('Federal', 'State', 'Local', 
                              'Interpersonal I', 'Interpersonal II'),
       title = 'Trust DV')

coefplot(type_vote1,main = '',
                   col = 1,
                   drop = c('Time', 'Black', 'Female', 'Class', 'Republican', 'PolInterest'),
                   sep = 0.5,
                   value.lab = 'Estimate with __ci__ Confidence Interval')
legend('bottomright', 
         col = 1, pch = c(20, 17),
       cex = 0.7,
       legend = c('2016 Voting', '2020 Voting'),
       title = 'Turnout DV')

coefplot(type_vote4,main = '',
                   col = 1,
                   drop = c('Time', 'Black', 'Female', 'Class', 'Republican', 'PolInterest'),
                   sep = 0.5,
                   value.lab = 'Estimate with __ci__ Confidence Interval')
legend('bottomright', 
         col = 1, pch = c(20, 17),
       cex = 0.7,
       legend = c('2016 Voting', '2020 Voting'),
       title = 'Turnout DV')

#### Analysis by Study

m_trst2 <- feols(c(trstgov, trststate, trstloc, trstgen, trstgen2) ~ PTS + PTG + Time +
                  Black + Female + Class + Republican + PolInterest | study,
                weights = ~weight,
                data = subset(datp, study == 'Two'))
summary(m_trst2)

m_trst3 <- feols(c(trstgov, trststate, trstloc, trstgen, trstgen2) ~ PTS + PTG + Time +
                  Black + Female + Class + Republican + PolInterest | study,
                weights = ~weight,
                data = subset(datp, study == 'Three'))
summary(m_trst3)

m_trst4 <- feols(c(trstgov, trststate, trstloc, trstgen, trstgen2) ~ PTS + PTG + Time +
                  Black + Female + Class + Republican + PolInterest | study,
                weights = ~weight,
                data = subset(datp, study == 'Four'))
summary(m_trst4)

m_trst5 <- feols(c(trstgov, trststate, trstloc, trstgen, trstgen2) ~ PTS + PTG + Time +
                  Black + Female + Class + Republican + PolInterest | study,
                weights = ~weight,
                data = subset(datp, study == 'Five'))
summary(m_trst5)


modelsummary(models = list('Federal'  = m_trst2[1],
                           'State'  = m_trst2[2],
                           'Local'  = m_trst2[3],
                           'Interpersonal I'  = m_trst2[4],
                           'Interpersonal II'  = m_trst2[5]),
             output = 'latex',
             stars = c('*' = .05),
             gof_omit = 'AIC|BIC|F|RMSE|Log.Lik.',
             escape = FALSE)

modelsummary(models = list('Federal'  = m_trst3[1],
                           'State'  = m_trst3[2],
                           'Local'  = m_trst3[3],
                           'Interpersonal I'  = m_trst3[4],
                           'Interpersonal II'  = m_trst3[5]),
             output = 'latex',
             stars = c('*' = .05),
             gof_omit = 'AIC|BIC|F|RMSE|Log.Lik.',
             escape = FALSE)

modelsummary(models = list('Federal'  = m_trst4[1],
                           'State'  = m_trst4[2],
                           'Local'  = m_trst4[3],
                           'Interpersonal I'  = m_trst4[4],
                           'Interpersonal II'  = m_trst4[5]),
             output = 'latex',
             stars = c('*' = .05),
             gof_omit = 'AIC|BIC|F|RMSE|Log.Lik.',
             escape = FALSE)

modelsummary(models = list('Federal'  = m_trst5[1],
                           'State'  = m_trst5[2],
                           'Local'  = m_trst5[3],
                           'Interpersonal I'  = m_trst5[4],
                           'Interpersonal II'  = m_trst5[5]),
             output = 'latex',
             stars = c('*' = .05),
             gof_omit = 'AIC|BIC|F|RMSE|Log.Lik.',
             escape = FALSE)


### Duration Effects

m_time1 <- feols(c(trstgov, trststate, trstloc, trstgen, trstgen2) ~ 
                  PTS + PTG +
                  Black + Female + Class + Republican + PolInterest | study,
                weights = ~weight,
                data = subset(datp, Time2 == 1))
summary(m_time1)

m_time2 <- feols(c(trstgov, trststate, trstloc, trstgen, trstgen2) ~ 
                  PTS + PTG +
                  Black + Female + Class + Republican + PolInterest | study,
                weights = ~weight,
                data = subset(datp, Time2 == 2))
summary(m_time2)

m_time3 <- feols(c(trstgov, trststate, trstloc, trstgen, trstgen2) ~ 
                  PTS + PTG +
                  Black + Female + Class + Republican + PolInterest | study,
                weights = ~weight,
                data = subset(datp, Time2 == 3))
summary(m_time3)

m_time4 <- feols(c(trstgov, trststate, trstloc, trstgen, trstgen2) ~ 
                  PTS + PTG +
                  Black + Female + Class + Republican + PolInterest | study,
                weights = ~weight,
                data = subset(datp, Time2 == 4))
summary(m_time4)

# modelsummary(models = list('Federal'  = m_trst5[1],
#                            'State'  = m_trst5[2],
#                            'Local'  = m_trst5[3],
#                            'Interpersonal I'  = m_trst5[4],
#                            'Interpersonal II'  = m_trst5[5]),
#              output = 'latex',
#              stars = c('*' = .05),
#              gof_omit = 'AIC|BIC|F|RMSE|Log.Lik.',
#              escape = FALSE)

coefplot(m_time1,main = '',
                   col = 1,
                   drop = c('Black', 'Female', 'Class', 'Republican', 'PolInterest'),
                   value.lab = 'Estimate with __ci__ Confidence Interval')
legend('bottomright', xpd = T,
         col = 1, pch = c(20, 17, 15, 21, 24),
       cex = 0.7,
       legend = c('Federal', 'State', 'Local', 
                              'Interpersonal I', 'Interpersonal II'),
       title = 'Trust DV')

coefplot(m_time2,main = '',
                   col = 1,
                   drop = c('Black', 'Female', 'Class', 'Republican', 'PolInterest'),
                   value.lab = 'Estimate with __ci__ Confidence Interval')
legend('bottomright', xpd = T,
         col = 1, pch = c(20, 17, 15, 21, 24),
       cex = 0.7,
       legend = c('Federal', 'State', 'Local', 
                              'Interpersonal I', 'Interpersonal II'),
       title = 'Trust DV')

coefplot(m_time3,main = '',
                   col = 1,
                   drop = c('Black', 'Female', 'Class', 'Republican', 'PolInterest'),
                   value.lab = 'Estimate with __ci__ Confidence Interval')
legend('bottomright', xpd = T,
         col = 1, pch = c(20, 17, 15, 21, 24),
       cex = 0.7,
       legend = c('Federal', 'State', 'Local', 
                              'Interpersonal I', 'Interpersonal II'),
       title = 'Trust DV')

coefplot(m_time4,main = '',
                   col = 1,
                   drop = c('Black', 'Female', 'Class', 'Republican', 'PolInterest'),
                   value.lab = 'Estimate with __ci__ Confidence Interval')
legend('bottomright', xpd = T,
         col = 1, pch = c(20, 17, 15, 21, 24),
       cex = 0.7,
       legend = c('Federal', 'State', 'Local', 
                              'Interpersonal I', 'Interpersonal II'),
       title = 'Trust DV')


#### Resource Deprivation 

m_resource <- feols(c(trstgov, trststate, trstloc, trstgen, trstgen2) ~ 
                  PTS*lower + PTG +
                  Time + Black + Female + Republican + PolInterest | study,
                weights = ~weight,
                data = datp)
summary(m_resource)

modelsummary(models = list('Federal'  = m_resource[1],
                           'State'  = m_resource[2],
                           'Local'  = m_resource[3],
                           'Interpersonal I'  = m_resource[4],
                           'Interpersonal II'  = m_resource[5]),
             output = 'latex',
             stars = c('*' = .05),
             gof_omit = 'AIC|BIC|F|RMSE|Log.Lik.',
             escape = FALSE)

# mediation analysis
library(mediation)

fit.med <- lm(trstgov ~ PTS,
              data = datp)
fit.dv <- lm(presvote20 ~ PTS + trstgov,
              data = datp)

med1 <- mediate(fit.med, fit.dv,
                treat = 'PTS',
                mediator = 'trstgov',
                boot = T)

fit.med2 <- lm(trstgen2 ~ PTS,
              data = datp)
fit.dv2 <- lm(presvote20 ~ PTS + trstgen2,
              data = datp)

med2 <- mediate(fit.med2, fit.dv2,
                treat = 'PTS',
                mediator = 'trstgen2',
                boot = T)

fit.med3 <- lm(trstgen2 ~ PTG,
              data = datp)
fit.dv3 <- lm(presvote20 ~ PTG + trstgen2,
              data = datp)

med3 <- mediate(fit.med3, fit.dv3,
                treat = 'PTG',
                mediator = 'trstgen2',
                boot = T)

fit.med4 <- lm(trstgen2 ~ PTG,
              data = datp)
fit.dv4 <- lm(presvote20 ~ PTG + trstgen2,
              data = datp)

med4 <- mediate(fit.med4, fit.dv4,
                treat = 'PTG',
                mediator = 'trstgen2',
                boot = T)

summary(med1); summary(med2); summary(med3); summary(med4)

modelsummary(models = list('PTS Federal'  = med1,
                           'PTG Interpersonal II' = med2,
                           'PTS Federal'  = med3,
                           'PTG Interpersonal II' = med4),
             output = 'latex',
             # stars = c('*' = .05),
             # gof_omit = 'AIC|BIC|F|RMSE|Log.Lik.',
             escape = FALSE)

modelsummary(med1)

plot(med1)

```

### SEM

```{r}

sem1 <-'
# factor construction
  ptgi1 =~ ptgi_1 + ptgi_2 + ptgi_3 + ptgi_4 + ptgi_5
  ptsi1 =~ ptsi_1 + ptsi_2 + ptsi_3 + ptsi_4 + ptsi_5
# models
  ptgi1 ~ black + hisp + female + lower + gop + newsint
  ptsi1 ~ black + hisp + female + lower + gop + newsint
  trstgov ~ ptgi1 + ptsi1 + black + hisp + female + lower + gop + newsint
  trststate ~ ptgi1 + ptsi1 + black + hisp + female + lower + gop + newsint
  trstloc ~ ptgi1 + ptsi1 + black + hisp + female + lower + gop + newsint
  trstgen ~ ptgi1 + ptsi1 + black + hisp + female + lower + gop + newsint
  trstgen2 ~ ptgi1 + ptsi1 + black + hisp + female + lower + gop + newsint
# residual correlations
  black ~~ female
  black ~~ lower
  female ~~ lower
  ptgi1 ~~ ptsi1
  trstgov ~~ trststate
  trstgov ~~ trstgen
  trststate ~~ trstgen
'

fit1 <- sem(sem1, missing = 'fiml', sampling.weights = 'weight', cluster = 'study',
            data = datp)
summary(fit1)

p_sem1 <- graph_sem(model = fit1)

## analyzing by type of traumatic event

fit2 <- sem(sem1, missing = 'fiml', sampling.weights = 'weight', cluster = 'study',
            group = 'Type',
            data = datp)
summary(fit2)

```

# Analysis

## SEM Approach

```{r}

sem1 <-'
# factor construction
  ptgi1 =~ ptgi_1 + ptgi_2 + ptgi_3 + ptgi_4 + ptgi_5
  ptsi1 =~ ptsi_1 + ptsi_2 + ptsi_3 + ptsi_4 + ptsi_5
# models
  ptgi ~ black + female + income
  ptsi ~ black + female + income
  trstgov ~ ptgi + ptsi + black + female + income + gop
  trststate ~ ptgi + ptsi + black + female + income + gop
  trstgen ~ ptgi + ptsi + black + female + income + gop
# residual correlations
  black ~~ female
  black ~~ income
  female ~~ income
  ptgi1 ~~ ptsi1
  trstgov ~~ trststate
  trstgov ~~ trstgen
  trststate ~~ trstgen
'

fit1 <- sem(sem1, missing = 'fiml', data = dat1)
summary(fit1)

sem2 <-'
# factor construction
  ptgi1 =~ ptgi_1 + ptgi_2 + ptgi_3 + ptgi_4 + ptgi_5
  ptsi1 =~ ptsi_1 + ptsi_2 + ptsi_3 + ptsi_4 + ptsi_5
# models
  ptgi ~ black + female + income
  ptsi ~ black + female + income
  trstgov ~ ptgi + ptsi + black + female + income + gop
  trststate ~ ptgi + ptsi + black + female + income + gop
  trstloc ~ ptgi + ptsi + black + female + income + gop
  trstgen ~ ptgi + ptsi + black + female + income + gop
  trstgen2 ~ ptgi + ptsi + black + female + income + gop
# residual correlations
  black ~~ female
  black ~~ income
  female ~~ income
  ptgi1 ~~ ptsi1
  trstgov ~~ trststate
  trstgov ~~ trstgen
  trststate ~~ trstgen
'

fit2 <- sem(sem2, missing = 'fiml', sampling.weights = 'weight', data = dat2)
summary(fit2)

sem3 <-'
# factor construction
  ptgi1 =~ ptgi_1 + ptgi_2 + ptgi_3 + ptgi_4 + ptgi_5
  ptsi1 =~ ptsi_1 + ptsi_2 + ptsi_3 + ptsi_4 + ptsi_5
# models
  ptgi ~ black + female + income
  ptsi ~ black + female + income
  trstgov ~ ptgi + ptsi + black + female + income + gop
  trststate ~ ptgi + ptsi + black + female + income + gop
  trstloc ~ ptgi + ptsi + black + female + income + gop
  trstgen ~ ptgi + ptsi + black + female + income + gop
  trstgen2 ~ ptgi + ptsi + black + female + income + gop
# residual correlations
  black ~~ female
  black ~~ income
  female ~~ income
  ptgi1 ~~ ptsi1
  trstgov ~~ trststate
  trstgov ~~ trstgen
  trststate ~~ trstgen
'

fit3 <- sem(sem3, missing = 'fiml', sampling.weights = 'weight', data = dat3)
summary(fit3)

sem4 <-'
# factor construction
  ptgi1 =~ ptgi_1 + ptgi_2 + ptgi_3 + ptgi_4 + ptgi_5
  ptsi1 =~ ptsi_1 + ptsi_2 + ptsi_3 + ptsi_4 + ptsi_5
# models
  ptgi ~ black + female + income
  ptsi ~ black + female + income
  trstgov ~ ptgi + ptsi + black + female + income + gop
  trststate ~ ptgi + ptsi + black + female + income + gop
  trstloc ~ ptgi + ptsi + black + female + income + gop
  trstgen ~ ptgi + ptsi + black + female + income + gop
  trstgen2 ~ ptgi + ptsi + black + female + income + gop
# residual correlations
  black ~~ female
  black ~~ income
  female ~~ income
  ptgi1 ~~ ptsi1
  trstgov ~~ trststate
  trstgov ~~ trstgen
  trststate ~~ trstgen
'

dat4$weight <- as.numeric(dat4$weight)
dat4 <- dat4 %>% drop_na(weight)
dat4 <- dat4 %>% drop_na(gop)


fit4 <- sem(sem4, missing = 'fiml', #sampling.weights = 'weight',
            data = dat4)
summary(fit4)

sem5 <-'
# factor construction
  ptgi1 =~ ptgi_1 + ptgi_2 + ptgi_3 + ptgi_4 + ptgi_5
  ptsi1 =~ ptsi_1 + ptsi_2 + ptsi_3 + ptsi_4 + ptsi_5
# models
  ptgi ~ black + female + income
  ptsi ~ black + female + income
  trstgov ~ ptgi + ptsi + black + female + income + gop
  trststate ~ ptgi + ptsi + black + female + income + gop
  trstloc ~ ptgi + ptsi + black + female + income + gop
  trstgen ~ ptgi + ptsi + black + female + income + gop
  trstgen2 ~ ptgi + ptsi + black + female + income + gop
# residual correlations
  black ~~ female
  black ~~ income
  female ~~ income
  ptgi1 ~~ ptsi1
  trstgov ~~ trststate
  trstgov ~~ trstgen
  trststate ~~ trstgen
'

fit5 <- sem(sem5, missing = 'fiml', #sampling.weights = 'weight',
            data = dat5)
summary(fit5)

race1 <-'
# factor construction
  ptgi1 =~ ptgi_1 + ptgi_2 + ptgi_3 + ptgi_4 + ptgi_5
  ptsi1 =~ ptsi_1 + ptsi_2 + ptsi_3 + ptsi_4 + ptsi_5
# models
  trstgov ~ ptgi + ptsi + gop + female + income
  trststate ~ ptgi + ptsi + gop + female + income
  trstgen ~ ptgi + ptsi + gop + female + income
# residual correlations
  female ~~ income
  gop ~~ income
  ptgi1 ~~ ptsi1
  trstgov ~~ trststate
  trstgov ~~ trstgen
  trststate ~~ trstgen
'

race2 <-'
# factor construction
  ptgi1 =~ ptgi_1 + ptgi_2 + ptgi_3 + ptgi_4 + ptgi_5
  ptsi1 =~ ptsi_1 + ptsi_2 + ptsi_3 + ptsi_4 + ptsi_5
# models
  trstgov ~ ptgi + ptsi + gop + female + income
  trststate ~ ptgi + ptsi + gop + female + income
  trstloc ~ ptgi + ptsi + gop + female + income
  trstgen ~ ptgi + ptsi + gop + female + income
  trstgen2 ~ ptgi + ptsi + gop + female + income
# residual correlations
  female ~~ income
  gop ~~ income
  ptgi1 ~~ ptsi1
  trstgov ~~ trststate
  trstgov ~~ trstgen
  trststate ~~ trstgen
'

# Constraining and Loosening across groups to test for differences

# study I
fit_r1.1 <- sem(race1, missing = 'fiml', group = 'black', #sampling.weights = 'weight',
            data = dat1)
fit_r1.2 <- sem(race1, missing = 'fiml', group = 'black', #sampling.weights = 'weight',
            data = dat1, group.equal = "loadings")
fit_r1.3 <- sem(race1, missing = 'fiml', group = 'black', #sampling.weights = 'weight',
            data = dat1, group.equal = c("intercepts", "loadings"))
lavTestLRT(fit_r1.1, fit_r1.2, fit_r1.3) # the groups are statistically significantly different

# pool all YouGov surveys 
dat_mg <- do.call("rbind", list(dat2, dat3, dat4, dat5))

# study II
fit_r1 <- sem(race2, missing = 'fiml', group = 'black', #sampling.weights = 'weight',
            data = dat_mg)
fit_r2 <- sem(race2, missing = 'fiml', group = 'black', #sampling.weights = 'weight',
            data = dat_mg, group.equal = "loadings")
fit_r3 <- sem(race2, missing = 'fiml', group = 'black', #sampling.weights = 'weight',
            data = dat_mg, group.equal = c("intercepts", "loadings"))
lavTestLRT(fit_r1, fit_r2, fit_r3) # the groups are statistically significantly different

summary(fit_r1)

modelsummary(models = list('Pooled Studies II-V' = fit_r1),
             output = 'latex',
             stars = c('*' = 0.05),
             gof_omit = 'AIC|BIC|F|Log.Lik.',
             escape = F)


## Gender

woman2 <-'
# factor construction
  ptgi1 =~ ptgi_1 + ptgi_2 + ptgi_3 + ptgi_4 + ptgi_5
  ptsi1 =~ ptsi_1 + ptsi_2 + ptsi_3 + ptsi_4 + ptsi_5
# models
  trstgov ~ ptgi + ptsi + gop + black + income
  trststate ~ ptgi + ptsi + gop + black + income
  trstloc ~ ptgi + ptsi + gop + black + income
  trstgen ~ ptgi + ptsi + gop + black + income
  trstgen2 ~ ptgi + ptsi + gop + black + income
# residual correlations
  black ~~ income
  gop ~~ income
  ptgi1 ~~ ptsi1
  trstgov ~~ trststate
  trstgov ~~ trstgen
  trststate ~~ trstgen
'

# study II
fit_w1 <- sem(woman2, missing = 'fiml', group = 'female', #sampling.weights = 'weight',
            data = dat_mg)
fit_w2 <- sem(woman2, missing = 'fiml', group = 'female', #sampling.weights = 'weight',
            data = dat_mg, group.equal = "loadings")
fit_w3 <- sem(woman2, missing = 'fiml', group = 'female', #sampling.weights = 'weight',
            data = dat_mg, group.equal = c("intercepts", "loadings"))
lavTestLRT(fit_w1, fit_w2, fit_w3) # the groups are statistically significantly different

summary(fit_w1)

modelsummary(models = list('Pooled Studies II-V' = fit_r1),
             output = 'latex',
             coef_omit = c(1:10, 18:20, 23:25, 28:30, 33:37, 39:88),
             # coef_rename = c('Non-Black' = 'Group1',
             #                  'Black' = 'Group 2'),
             stars = c('*' = 0.05),
             gof_omit = 'AIC|BIC|F|Log.Lik.',
             escape = F)

```

## Modeling Voting

### Path Model Approach

Stress consistently associated with increased trust in government and decreased likelihood of voting in 2016 across samples 2-5. Significant PTSI--2016 voting in samples 2, 3, 5 and PTSI--2020 voting in samples 2, 3, 4, 5. PTSI--trust in federal government is significant in 3, 4, 5. PTGI is not significantly associated with trust in the federal government or with voting in 2016 or 2020 in any samples.

Multiple group SEM to identify if effect is divergent for Black and White respondents. Due to representativeness of sample, number of Black respondents is low. I pool samples 2 and 3 as well as 4 and 5 given they were fielded at the same time and to increase sample size of Black respondents. No significant differences between Black and White respondents in terms of these associations. 

I then check for female versus male respondents and income groups. 

For income levels, I bisect the samples comparing those who make less than \$50,000 per year in family income to those who make \$50,000 or more per year in family income. In terms of significant associations: PTSI--2016 voting

```{r Path Model--Voting}

## PTSI

med_1s <- '# direct effect
        presvote16 ~ c*ptsi
        # mediator
        trstgov ~ a*ptsi
        presvote16 ~ b*trstgov
        # indirect effect
        ab := a*b
        # total effect
        total := c + (a*b)
        '
med_2s <- '# direct effect
        presvote20 ~ c*ptsi
        # mediator
        trstgov ~ a*ptsi
        presvote20 ~ b*trstgov
        # indirect effect
        ab := a*b
        # total effect
        total := c + (a*b)
        '

med2_1s <- sem(med_1s, missing = 'fiml', sampling.weights = 'weight', data = dat2)

med2_2s <- sem(med_2s, missing = 'fiml', sampling.weights = 'weight', data = dat2)

med3_1s <- sem(med_1s, missing = 'fiml', sampling.weights = 'weight', data = dat3)

med3_2s <- sem(med_2s, missing = 'fiml', sampling.weights = 'weight', data = dat3)

med4_1s <- sem(med_1s, missing = 'fiml', sampling.weights = 'weight', data = dat4)

med4_2s <- sem(med_2s, missing = 'fiml', sampling.weights = 'weight', data = dat4)

med5_1s <- sem(med_1s, missing = 'fiml', sampling.weights = 'weight', data = dat5)

med5_2s <- sem(med_2s, missing = 'fiml', sampling.weights = 'weight', data = dat5)

summary(med2_1s); summary(med2_2s); summary(med3_1s); summary(med3_2s); summary(med4_1s); summary(med4_2s); summary(med5_1s); summary(med5_2s)

modelsummary(models = list('2016'  = med2_1s,
                           '2020'  = med2_2s,
                           '2016'  = med3_1s,
                           '2020'  = med3_2s,
                           '2016'  = med4_1s,
                           '2020'  = med4_2s,
                           '2016'  = med5_1s,
                           '2020'  = med5_2s),
             output = 'latex',
             stars = c('*' = .05),
             gof_omit = 'AIC|BIC|F|RMSE|Log.Lik.',
             escape = FALSE)

## PTGI

med_1g <- '# direct effect
        presvote16 ~ c*ptgi
        # mediator
        trstgov ~ a*ptgi
        presvote16 ~ b*trstgov
        # indirect effect
        ab := a*b
        # total effect
        total := c + (a*b)
        '
med_2g <- '# direct effect
        presvote20 ~ c*ptgi
        # mediator
        trstgov ~ a*ptgi
        presvote20 ~ b*trstgov
        # indirect effect
        ab := a*b
        # total effect
        total := c + (a*b)
        '

med2_1g <- sem(med_1g, missing = 'fiml', sampling.weights = 'weight', data = dat2)

med2_2g <- sem(med_2g, missing = 'fiml', sampling.weights = 'weight', data = dat2)

med3_1g <- sem(med_1g, missing = 'fiml', sampling.weights = 'weight', data = dat3)

med3_2g <- sem(med_2g, missing = 'fiml', sampling.weights = 'weight', data = dat3)

med4_1g <- sem(med_1g, missing = 'fiml', sampling.weights = 'weight', data = dat4)

med4_2g <- sem(med_2g, missing = 'fiml', sampling.weights = 'weight', data = dat4)

med5_1g <- sem(med_1g, missing = 'fiml', sampling.weights = 'weight', data = dat5)

med5_2g <- sem(med_2g, missing = 'fiml', sampling.weights = 'weight', data = dat5)

summary(med2_1g); summary(med2_2g); summary(med3_1g); summary(med3_2g); summary(med4_1g); summary(med4_2g); summary(med5_1g); summary(med5_2g)

modelsummary(models = list('2016'  = med2_1g,
                           '2020'  = med2_2g,
                           '2016'  = med3_1g,
                           '2020'  = med3_2g,
                           '2016'  = med4_1g,
                           '2020'  = med4_2g,
                           '2016'  = med5_1g,
                           '2020'  = med5_2g),
             output = 'latex',
             stars = c('*' = .05),
             gof_omit = 'AIC|BIC|F|RMSE|Log.Lik.',
             escape = FALSE)

# ### Multiple Group SEM
# 
# #### Race
# 
# dat23 <- rbind(dat2, dat3)
# dat45 <- rbind(dat4, dat5)
# 
# grp2_1s <- sem(med_1s, missing = 'fiml', sampling.weights = 'weight', group = 'black', data = dat23)
# 
# grp2_2s <- sem(med_2s, missing = 'fiml', sampling.weights = 'weight', group = 'black', data = dat23)
# 
# grp4_1s <- sem(med_1s, missing = 'fiml', sampling.weights = 'weight', group = 'black', data = dat45)
# 
# grp4_2s <- sem(med_2s, missing = 'fiml', sampling.weights = 'weight', group = 'black', data = dat45)
# 
# summary(grp2_1s); summary(grp2_2s); summary(grp4_1s); summary(grp4_2s)
# 
# grp2_1g <- sem(med_1g, missing = 'fiml', sampling.weights = 'weight', group = 'black', data = dat23)
# 
# grp2_2g <- sem(med_2g, missing = 'fiml', sampling.weights = 'weight', group = 'black', data = dat23)
# 
# grp4_1g <- sem(med_1g, missing = 'fiml', sampling.weights = 'weight', group = 'black', data = dat45)
# 
# grp4_2g <- sem(med_2g, missing = 'fiml', sampling.weights = 'weight', group = 'black', data = dat45)
# 
# summary(grp2_1g); summary(grp2_2g); summary(grp4_1g); summary(grp4_2g)

# #### Woman
# 
# dat23 <- rbind(dat2, dat3)
# dat45 <- rbind(dat4, dat5)
# 
# grp2_1s <- sem(med_1s, missing = 'fiml', sampling.weights = 'weight', group = 'black', data = dat23)
# 
# grp2_2s <- sem(med_2s, missing = 'fiml', sampling.weights = 'weight', group = 'black', data = dat23)
# 
# grp4_1s <- sem(med_1s, missing = 'fiml', sampling.weights = 'weight', group = 'black', data = dat45)
# 
# grp4_2s <- sem(med_2s, missing = 'fiml', sampling.weights = 'weight', group = 'black', data = dat45)
# 
# summary(grp2_1s); summary(grp2_2s); summary(grp4_1s); summary(grp4_2s)
# 
# grp2_1g <- sem(med_1g, missing = 'fiml', sampling.weights = 'weight', group = 'black', data = dat23)
# 
# grp2_2g <- sem(med_2g, missing = 'fiml', sampling.weights = 'weight', group = 'black', data = dat23)
# 
# grp4_1g <- sem(med_1g, missing = 'fiml', sampling.weights = 'weight', group = 'black', data = dat45)
# 
# grp4_2g <- sem(med_2g, missing = 'fiml', sampling.weights = 'weight', group = 'black', data = dat45)
# 
# summary(grp2_1g); summary(grp2_2g); summary(grp4_1g); summary(grp4_2g)
# 
# #### Income
# 
# dat2$class <- ifelse(dat2$income < 0.3, 0, 1)
# dat3$class <- ifelse(dat3$income < 0.3, 0, 1)
# dat4$class <- ifelse(dat4$income < 0.3, 0, 1)
# dat5$class <- ifelse(dat5$income < 0.3, 0, 1)
# 
# grp2_1s <- sem(med_1s, missing = 'fiml', sampling.weights = 'weight', group = 'class', data = dat2)
# 
# grp2_2s <- sem(med_2s, missing = 'fiml', sampling.weights = 'weight', group = 'class', data = dat2)
# 
# grp3_1s <- sem(med_1s, missing = 'fiml', sampling.weights = 'weight', group = 'class', data = dat3)
# 
# grp3_2s <- sem(med_2s, missing = 'fiml', sampling.weights = 'weight', group = 'class', data = dat3)
# 
# grp4_1s <- sem(med_1s, missing = 'fiml', sampling.weights = 'weight', group = 'class', data = dat4)
# 
# grp4_2s <- sem(med_2s, missing = 'fiml', sampling.weights = 'weight', group = 'class', data = dat4)
# 
# grp5_1s <- sem(med_1s, missing = 'fiml', sampling.weights = 'weight', group = 'class', data = dat5)
# 
# grp5_2s <- sem(med_2s, missing = 'fiml', sampling.weights = 'weight', group = 'class', data = dat5)
# 
# summary(grp2_1s); summary(grp2_2s); summary(grp3_1s); summary(grp3_2s); summary(grp4_1s); summary(grp4_2s); summary(grp5_1s); summary(grp5_2s)
# 
# grp2_1g <- sem(med_1g, missing = 'fiml', sampling.weights = 'weight', group = 'black', data = dat23)
# 
# grp2_2g <- sem(med_2g, missing = 'fiml', sampling.weights = 'weight', group = 'black', data = dat23)
# 
# grp4_1g <- sem(med_1g, missing = 'fiml', sampling.weights = 'weight', group = 'black', data = dat45)
# 
# grp4_2g <- sem(med_2g, missing = 'fiml', sampling.weights = 'weight', group = 'black', data = dat45)
# 
# summary(grp2_1g); summary(grp2_2g); summary(grp4_1g); summary(grp4_2g)



```

## Temporal Proximity

```{r Temporal Proximity}

# Trust

#self exposure time
t2_t <- glm(trstgov ~ tr_exp_t + age + edu + female, data = dat2, family = 'binomial', weights = weight)
t3_t <- glm(trstgov ~ tr_exp_t + age + edu + female, data = dat3, family = 'binomial', weights = weight)
t4_t <- glm(trstgov ~ tr_exp_t + age + edu + female, data = dat4, family = 'binomial', weights = as.numeric(weight))
t5_t <- glm(trstgov ~ tr_exp_t + age + edu + female, data = dat5, family = 'binomial', weights = weight)

summary(t2_t); summary(t3_t); summary(t4_t); summary(t5_t)

# close friend/family member exposure time
t2_rt <- glm(trstgov ~ rtr_exp_t + age + edu + female, data = dat2, family = 'binomial', weights = weight)
t3_rt <- glm(trstgov ~ rtr_exp_t + age + edu + female, data = dat3, family = 'binomial', weights = weight)
t4_rt <- glm(trstgov ~ rtr_exp_t + age + edu + female, data = dat4, family = 'binomial', weights = as.numeric(weight))
t5_rt <- glm(trstgov ~ rtr_exp_t + age + edu + female, data = dat5, family = 'binomial', weights = weight)

summary(t2_rt); summary(t3_rt); summary(t4_rt); summary(t5_rt)

# turnout
# trauma exposure time
t2_tr <- glm(presvote20 ~ tr_exp_t + age + edu + female, data = dat2, family = 'binomial', weights = weight)
t3_tr <- glm(presvote20 ~ tr_exp_t + age + edu + female, data = dat3, family = 'binomial', weights = weight)
t4_tr <- glm(presvote20 ~ tr_exp_t + age + edu + female, data = dat4, family = 'binomial', weights = as.numeric(weight))
t5_tr <- glm(presvote20 ~ tr_exp_t + age + edu + female, data = dat5, family = 'binomial', weights = weight)

summary(t2_tr); summary(t3_tr); summary(t4_tr); summary(t5_tr)

# trauma exposure time interacted with PTSI
t2_tr2 <- glm(presvote20 ~ tr_exp_t*ptsi + age + edu + female, data = dat2, family = 'binomial', weights = weight)
t3_tr2 <- glm(presvote20 ~ tr_exp_t*ptsi + age + edu + female, data = dat3, family = 'binomial', weights = weight)
t4_tr2 <- glm(presvote20 ~ tr_exp_t*ptsi + age + edu + female, data = dat4, family = 'binomial', weights = as.numeric(weight))
t5_tr2 <- glm(presvote20 ~ tr_exp_t*ptsi + age + edu + female, data = dat5, family = 'binomial', weights = weight)

summary(t2_tr2); summary(t3_tr2); summary(t4_tr2); summary(t5_tr2)

# trauma exposure time interacted with PTGI
t2_tr3 <- glm(presvote20 ~ tr_exp_t*ptgi + age + edu + female, data = dat2, family = 'binomial', weights = weight)
t3_tr3 <- glm(presvote20 ~ tr_exp_t*ptgi + age + edu + female, data = dat3, family = 'binomial', weights = weight)
t4_tr3 <- glm(presvote20 ~ tr_exp_t*ptgi + age + edu + female, data = dat4, family = 'binomial', weights = as.numeric(weight))
t5_tr3 <- glm(presvote20 ~ tr_exp_t*ptgi + age + edu + female, data = dat5, family = 'binomial', weights = weight)

summary(t2_tr3); summary(t3_tr3); summary(t4_tr3); summary(t5_tr3)

# close friend/family member trauma time
t2_rtr <- glm(presvote20 ~ rtr_exp_t + age + edu + female, data = dat2, family = 'binomial', weights = weight)
t3_rtr <- glm(presvote20 ~ rtr_exp_t + age + edu + female, data = dat3, family = 'binomial', weights = weight)
t4_rtr <- glm(presvote20 ~ rtr_exp_t + age + edu + female, data = dat4, family = 'binomial', weights = as.numeric(weight))
t5_rtr <- glm(presvote20 ~ rtr_exp_t + age + edu + female, data = dat5, family = 'binomial', weights = weight)

summary(t2_rtr); summary(t3_rtr); summary(t4_rtr); summary(t5_rtr)

```


### Regression Approach (OLS and Logit GLM)

```{r Regression--Voting}

m1_2 <- glm(presvote16 ~ ptgi*ptsi*trstgov, data = dat2, family = 'binomial', weights = weight)
summary(m1_2)

m2_2 <- glm(presvote20 ~ ptgi*ptsi*trstgov + presvote16, data = dat2, family = 'binomial', weights = weight)
summary(m2_2)

m1_3 <- glm(presvote16 ~ ptgi*ptsi*trstgov, data = dat3, family = 'binomial', weights = weight)
summary(m1_3)

m2_3 <- glm(presvote20 ~ ptgi*ptsi*trstgov + presvote16, data = dat3, family = 'binomial', weights = weight)
summary(m2_3)

m1_4 <- glm(presvote16 ~ ptgi*ptsi*trstgov, data = dat4, family = 'binomial', weights = weight)
summary(m1_4)

m2_4 <- glm(presvote20 ~ ptgi*ptsi*trstgov + presvote16, data = dat4, family = 'binomial', weights = weight)
summary(m2_4)

m1_5 <- glm(presvote16 ~ ptgi*ptsi*trstgov, data = dat5, family = 'binomial', weights = weight)
summary(m1_5)

m2_5 <- glm(presvote20 ~ ptgi*ptsi*trstgov + presvote16, data = dat5, family = 'binomial', weights = weight)
summary(m2_5)

# Trust in Federal Government, by party

m1_1 <- lm(trstgov ~ ptsi + ptgi + female + black + income + gop,
           data = dat1)

m1r_1 <- lm(trstgov ~ ptsi + ptgi + female + black + income,
           data = subset(dat1, gop == 1))
m1d_1 <- lm(trstgov ~ ptsi + ptgi + female + black + income,
           data = subset(dat1, gop == 0))
m1r_2 <- lm(trstgov ~ ptsi + ptgi + female + black + income,
           data = subset(dat2, gop == 1))
m1d_2 <- lm(trstgov ~ ptsi + ptgi + female + black + income,
           data = subset(dat2, gop == 0))
m1r_3 <- lm(trstgov ~ ptsi + ptgi + female + black + income,
           data = subset(dat3, gop == 1))
m1d_3 <- lm(trstgov ~ ptsi + ptgi + female + black + income,
           data = subset(dat3, gop == 0))

m2_1 <- lm(trststate ~ ptsi + ptgi + female + black + income + gop,
           data = dat1)

m3_1 <- glm(trstgen ~ ptsi + ptgi + female + black + income + gop,
           family = 'binomial', data = dat1)

```


